Introduccion

El modelamiento del comportamiento de diferentes variables es un tema que ha sido estudiado en sectores energéticos, industriales, económicos y financieros. De allí se comienza a apreciar tanto la importancia que tienen los datos hoy en día al igual que las técnicas utilizadas para su modelamiento. La estadística es una disciplina que se preocupa por la recolección, organización, interpretación y análisis de datos, que, según su aplicación puede traer un gran impacto en la industria al momento de la toma de decisiones. En particular, diferentes técnicas estadísticas han sido utilizadas en el sector financiero, las cuales permiten modelar comportamiento de los clientes, acciones, entre otras variables.

Diferentes industrias dentro de su funcionamiento, deben presentar ante la superintendencia información relacionando los gastos y ventas que presentaron anualmente. Además, se presume que del mercado colombiano, como en los procesos de las industria está esa interacción con todo el mercado, es posible pensar que exista una relación entre las diferentes variables macroeconomicas (ej. PIB, TRM, Balance Fiscal, Indice de Desempleo, etc.) y estos montos de costos y gastos de las empresas. Debido a la cantidad de información con la que se cuenta (información de costos y gastos para diferentes empresas en colombia entre los años 2016 y 2018), se sabe que no se cuenta con información suficiente para la construcción de un modelo por empresa que permita ver la relación existente entre las variables macroeconómicas y las variables asociadas a costos y gastos de ventas. Por lo anterior, es posible considerar un conjunto de datos como la consolidación de la información de todas las empresas junto con la información macroeconomica para los años en estudio, así buscando construir un modelo general para realizar la modelación de estas variables reportadas ante la superintendencia para un conjunto de empresas cuya industria sea similar.

Por lo tanto, para el conjunto de datos meniconado anteriormente, se buscará modelar la información de costos y gastos de venta a partir de las variables macroeconomicas disponibles, al igual que analizar si al realizar alguna transformación a dichas variables resulta relevante al momento de la creación del modelo. Se utilizarán modelos lineales, comenzando con la evaluación del modelo lineal general, hasta la aplicación de modelos de efectos mixtos. Esta última estructura de modelos, es bastante usado al momento de tener individuos que comparten la misma información pero tienen una salida diferente (en nuestro caso, todas las empresas comparten la misma información de las variables macroeconómicas), por lo que utilizar este tipo de modelos resulta de gran interés ya que permite modelar tanto efectos a nivel de individuo como agregando un efecto aleatorio.

Clasificación Industrial Seleccionada

El sector de la construcción es uno de los más relevantes en la economía colombiana. En nuestro contexto nacional el sector es considerado como uno de los más vitales para el desarrollo del país y representa uno de los más importantes rubros en materia de produccción interna componiendo cada año de 6 a 7 por ciento del producto interno bruto total y hasta un 7.1% del total de ocupados a nivel nacional. Dicho sector es caracterizado por sus fluctuaciones estacionarias fuertemente influenciadas por los planes de infraestructura de gran escala y los planes de gobierno. Respecto al último trimestre de 2019 tuvo un aumento de 3.4%, uno de los más altos a nivel de america latina.

Si bien las estimaciones globales para el año 2020 en Colombia para esta industria eran positivas, la coyuntura del COVID-19 ha de perturbar fuertemente el sector, afectando con alto impacto a los importadores de materiales y a la demanada frente a los retrasos en la ejecución de obras. Desde enero el sector de la producción de concreto ya estaba presentando caídas significativas de hasta un 8.3% respecto al año pasado en el mismo mes, por lo que se esperan peores resultados al cierre del segundo trimestre del año presente. Muchas compañías planearon incrementos en sus precios, con la esperanza de generar mayores ingresos, pero dichos planes han de ser postergados bajo la actual coyuntura. La demanda, el driver más relevante en la industria, claramente se ve desplazado efecto del impedimento de comercialización y la paralisis en el país bajo las medidas de cuarentena nacional. Un punto importante es que algunas compañías del sector podrán seguir con sus actividades producto de la inclusión de actividades de infraestructura como vital durante la crisis del COVID-19.

Frente a la incertidumbre que generan estos escenarios y como la dinámica particular de cada compañía que aporta al crecimiento de la industria total evoluciona en el tiempo se torna relevante construir análisis y modelos estadísticos robustos que respondan a las perturbaciones en el mercado y permitan obtener información valiosa para la toma de decisiones.

Consolidación de información

all_data = read.csv("PreparacionDatos/Datos_completos.csv", encoding="UTF-8")
NIT = all_data['NIT'] 
data = all_data[,-c(1,2)]

Análisis descriptivo

Damos una visualizacion inicial al conjunto de datos con el que se va a trabajar:

kable(head(round(data,2))) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "300px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
60626083 13978834 -0.10 0.14 0.11 0.07 0.09 0.06 0.52 0.35 0.27
3785517 60608 0.07 0.33 0.11 0.07 0.09 0.06 0.52 0.35 0.27
39634059 5645885 0.34 0.69 0.11 0.07 0.09 0.06 0.52 0.35 0.27
24304534 8977594 -0.05 -0.04 0.11 0.07 0.09 0.06 0.52 0.35 0.27
31280967 6702215 -0.05 -0.18 0.11 0.07 0.09 0.06 0.52 0.35 0.27
11886427 1050755 0.13 0.03 0.11 0.07 0.09 0.06 0.52 0.35 0.27

En la siguiente tabla, vemos información descriptiva de las variables que se considerarán en la etapa de modelamiento del trabajo, vemos así que contamos con un total de 11 variables, cada una con magnitudes diferentes. Vemos que algunos indicadores tienen valores entre 0 y 1 tal como lo es el desempleo, mientras que otras variables representan dinero, tal como la trm, el monto de gastos y ventas. Así, tenemos una idea inicial de las características del conjunto de datos, por lo que se aplicarán lo métodos correspondientes cuando sea necesario, si lo que se realizará es sensible a la escala de los datos.

summary_df = stat.desc(data)
kable(summary_df) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "300px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
nbr.val 9.900000e+01 9.900000e+01 99.0000000 99.0000000 99.0000000 99.0000000 99.0000000 99.0000000 99.0000000 99.0000000 99.0000000
nbr.null 0.000000e+00 0.000000e+00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
nbr.na 0.000000e+00 0.000000e+00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
min 5.913660e+05 3.436800e+04 -0.9289315 -0.7898731 -0.0334936 0.0656288 0.0922332 0.0318000 -0.2895537 -0.2740147 -0.2287460
max 2.119141e+08 3.541196e+07 1.0883866 6.1607837 0.1117616 0.0734318 0.0968134 0.0575000 0.5192874 0.3516454 0.2726036
range 2.113227e+08 3.537760e+07 2.0173181 6.9506568 0.1452552 0.0078030 0.0045802 0.0257000 0.8088411 0.6256601 0.5013496
sum 4.755036e+09 6.487101e+08 -3.7634585 6.2073859 2.6432290 6.9358217 9.3332892 4.2966000 3.0437504 7.4840537 -0.5776688
median 3.527843e+07 5.204078e+06 -0.0254994 0.0107538 0.0018298 0.0711158 0.0937804 0.0409000 -0.1374988 0.1491588 -0.0613627
mean 4.803067e+07 6.552627e+06 -0.0380147 0.0627009 0.0266993 0.0700588 0.0942756 0.0434000 0.0307450 0.0755965 -0.0058350
SE.mean 4.567598e+06 6.609484e+05 0.0279735 0.0665708 0.0062481 0.0003305 0.0001922 0.0010748 0.0354548 0.0263314 0.0210523
CI.mean.0.95 9.064249e+06 1.311631e+06 0.0555124 0.1321077 0.0123991 0.0006559 0.0003814 0.0021329 0.0703589 0.0522538 0.0417776
var 2.065432e+15 4.324843e+13 0.0774690 0.4387361 0.0038648 0.0000108 0.0000037 0.0001144 0.1244474 0.0686408 0.0438768
std.dev 4.544702e+07 6.576354e+06 0.2783325 0.6623716 0.0621675 0.0032887 0.0019121 0.0106940 0.3527710 0.2619939 0.2094678
coef.var 9.462084e-01 1.003621e+00 -7.3217010 10.5639935 2.3284333 0.0469421 0.0202817 0.2464057 11.4741097 3.4656883 -35.8982759

Además de la información descriptiva presentados en la tabla anterior, podemos ver para cada una de las variables, su distribución de forma visual con la ayuda de la creación de histogramas de frecuencia.

Vamos a obtener la informacion relacionada a las medidas de centralidad y dispersion del conjunto de datos. Inicialmente, presentamos información del vector de medias y medianas que describen la centralidad del conjunto de datos. Posteriormente, consideramos las matrices de Covarianza y Correlación para tener una intuición de la variabilidad de la información que consideramos.

Vector de Medias y Medianas

Media y Mediana
x
Costo.de.ventas 4.803067e+07
Gastos.de.ventas 6.552627e+06
Costo.de.ventas_dif -3.801470e-02
Gastos.de.ventas_dif 6.270090e-02
TRM 2.669930e-02
PIB 7.005880e-02
Desempleo 9.427560e-02
Inflacion 4.340000e-02
Tasa_Intervencion 3.074500e-02
Balance_CC 7.559650e-02
Balance_Fiscal -5.835000e-03
x
Costo.de.ventas 3.527843e+07
Gastos.de.ventas 5.204078e+06
Costo.de.ventas_dif -2.549940e-02
Gastos.de.ventas_dif 1.075380e-02
TRM 1.829800e-03
PIB 7.111580e-02
Desempleo 9.378040e-02
Inflacion 4.090000e-02
Tasa_Intervencion -1.374988e-01
Balance_CC 1.491588e-01
Balance_Fiscal -6.136270e-02

Análisis desde variables absolutas

Boxplot

par(mfrow=c(1,2))
bxplot_costos = boxplot(Costo.de.ventas~Year, data = all_data)
bxplot_gastos = boxplot(Gastos.de.ventas~Year, data = all_data)

### Outliers para costos
nits_tab = t(all_data[(all_data$Costo.de.ventas %in% bxplot_costos$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 860009694 900378893 860009694 900378893 860009694 900378893 
### Outliers para gastos
nits_tab = t(all_data[(all_data$Gastos.de.ventas %in% bxplot_gastos$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 801002644 801002644 800015615 801002644

Análisis desde variables diferenciales

Boxplot

par(mfrow=c(1,2))
bxplot_costos_dif = boxplot(Costo.de.ventas_dif~Year, data = all_data)
bxplot_gastos_dif = boxplot(Gastos.de.ventas_dif~Year, data = all_data)

### Outliers para costos dif

nits_tab = t(all_data[(all_data$Costo.de.ventas_dif %in% bxplot_costos_dif$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 800236890 890909034 900437650 800236890 890909034 900389088 800236890 
### Outliers para gastos dif

nits_tab = t(all_data[(all_data$Gastos.de.ventas_dif %in% bxplot_gastos_dif$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 800081030 890904459 900389088 800045720 800236890 890311366 890909034 900234565 900364670 900378893

Análisis:

Al observar la distribución del comportamiento de las variables de ventas y gastos para los 3 años en consideración, podemos ver que no hay mucha variación entre el comportamiento de estas a lo largo del tiempo, vemos que en particular, tanto el monto de costo y gasto promedio (mediana) es constante para los 3 años, aunque se nota un aumento en los costos en el año 2018 para el caso de algunas empresas. Ademas, vemos para que caso de los Costos, las mismas 2 empresas (relacionadas a los NITs 860009694 y 900378893) fueron las que presentaron un maypr valor de costos en los 3 periodos, por lo que el boxplot los considera registros outliers. Mirando la distribución de los gatos, notamos que una empresa tuvo los gastos mas altos en los 3 periodos (relacionadas al NIT 801002644) pero adicional, en el último año, la empresa con el NIT 800015615 tuvo un aumento en sus gastos, por lo que es considerado también como un outlier por el boxplot.

Pasando a la distribución de las variables de diferencia relativa de costos y gastos, vemos para ambos casos los registros outliers son aquellos que tienen un valor negativo en su diferencia o un incremento porcentual muy alto. Podemos ver que para el caso del 2016 en la variable costos, las empresas outliers fueron las correspondiente a los NITs 800236890, 890909034, 900437650, de estas, la única que no presenta una diferencia negativa fue la primera (800236890), por lo que puede indicar un aumento considerable (raro) en el mercado en término de los costos de venta (practicamente el doble del año anterior), mientras que todas las otras empresas presentaron una disminución considerable en este mismo concepto. Para el caso del 2017, tres empresas fueron consideradas como outlier, entre ellas, dos que se consideraron en el periodo anterior (800236890, 890909034), las cuales ambas reportaron una reducción de aproximadamente 65% en sus costos respecto al año anterior, mientras que la otra empresa considerada como ergistro raro (900389088) tuvo un aumento en sus costos de un 88% (casi el doble de los costos del año anterior). Para el último periodo, solo se encuentra una de las empresas mencionadas en los dos periodos anteriores (800236890), la cual vuelve a reportar una disminución considerable en sus costos.

Ahora mirando la variable de diferencia de gastos, solamente 2 empresas tienen comportamientos alto, presentando un aumento de casi el doble de los gastos anuales, estas empresas corresponden a 800081030 y 890904459 para los años 2016 y 2017 respectivamente, en cambio para el 2018, se tuvo variabilidad en los gastos de ventas de 8 empresas que son considerados anormales tanto por sus aumentos como disminuciones, auqellas empresas son 900389088, 800045720, 800236890, 890311366, 890909034, 900234565, 900364670 y 900378893. En aquellas empresas que tuvieron aumentos considerables, reportan casi hasta un cuarto del aumento de gastos en comparación al año anterior, mientras que aquellas que reportan una disminución, vemos que sus gastos se disminuyeron hasta el doble del periodo pasado.

Asi, estos cambios en las industrias en sus costos y gastos de ventas, pueden estar directamente relacionados con algunos eventos que hayan hecho la gran variación en la materia prima requerida para sus productos (tipos de mezlcas de cemento por ejemplo), además que el sector seleccionado depende bastante de los planes de infraestructura que se tengan, por lo que también es factible que en los periodos analizados, hayan ocurrido cambios en estos planes dentro de la contratación de cada empresa, y sean estas las causas que nos lleven a ver diferencias tan grandes de periodo a periodo (tanto positivas como negativas)

Matriz de Covarianzas

cov_data = cov(data)
kable(formatC(cov_data,format = "e", digits = 2)) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "480px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
Costo.de.ventas 2.07e+15 1.35e+14 2.34e+06 -1.56e+06 8.13e+04 5.00e+03 -2.64e+02 6.97e+03 3.12e+05 2.99e+04 2.89e+05
Gastos.de.ventas 1.35e+14 4.32e+13 3.58e+05 1.04e+05 3.18e+03 3.62e+02 1.48e+02 -3.37e+02 -2.24e+03 -2.08e+04 1.34e+04
Costo.de.ventas_dif 2.34e+06 3.58e+05 7.75e-02 -3.65e-02 1.14e-03 9.75e-05 2.24e-05 -3.04e-06 1.99e-03 -3.20e-03 4.39e-03
Gastos.de.ventas_dif -1.56e+06 1.04e+05 -3.65e-02 4.39e-01 8.23e-04 2.07e-04 1.46e-04 -5.03e-04 -1.05e-02 -2.03e-02 4.91e-03
TRM 8.13e+04 3.18e+03 1.14e-03 8.23e-04 3.86e-03 1.78e-04 -6.97e-05 5.52e-04 2.01e-02 9.34e-03 1.30e-02
PIB 5.00e+03 3.62e+02 9.75e-05 2.07e-04 1.78e-04 1.08e-05 -6.78e-07 1.56e-05 6.92e-04 7.94e-05 6.29e-04
Desempleo -2.64e+02 1.48e+02 2.24e-05 1.46e-04 -6.97e-05 -6.78e-07 3.66e-06 -1.92e-05 -5.82e-04 -5.01e-04 -2.02e-04
Inflacion 6.97e+03 -3.37e+02 -3.04e-06 -5.03e-04 5.52e-04 1.56e-05 -1.92e-05 1.14e-04 3.71e-03 2.61e-03 1.73e-03
Tasa_Intervencion 3.12e+05 -2.24e+03 1.99e-03 -1.05e-02 2.01e-02 6.92e-04 -5.82e-04 3.71e-03 1.24e-01 7.90e-02 6.44e-02
Balance_CC 2.99e+04 -2.08e+04 -3.20e-03 -2.03e-02 9.34e-03 7.94e-05 -5.01e-04 2.61e-03 7.90e-02 6.86e-02 2.69e-02
Balance_Fiscal 2.89e+05 1.34e+04 4.39e-03 4.91e-03 1.30e-02 6.29e-04 -2.02e-04 1.73e-03 6.44e-02 2.69e-02 4.39e-02

Matriz de Correlación

cor_data = cor(data)
kable(round(cor_data,2)) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "480px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
Costo.de.ventas 1.00 0.45 0.18 -0.05 0.03 0.03 0.00 0.01 0.02 0.00 0.03
Gastos.de.ventas 0.45 1.00 0.20 0.02 0.01 0.02 0.01 0.00 0.00 -0.01 0.01
Costo.de.ventas_dif 0.18 0.20 1.00 -0.20 0.07 0.11 0.04 0.00 0.02 -0.04 0.08
Gastos.de.ventas_dif -0.05 0.02 -0.20 1.00 0.02 0.10 0.12 -0.07 -0.04 -0.12 0.04
TRM 0.03 0.01 0.07 0.02 1.00 0.87 -0.59 0.83 0.92 0.57 1.00
PIB 0.03 0.02 0.11 0.10 0.87 1.00 -0.11 0.44 0.60 0.09 0.91
Desempleo 0.00 0.01 0.04 0.12 -0.59 -0.11 1.00 -0.94 -0.86 -1.00 -0.50
Inflacion 0.01 0.00 0.00 -0.07 0.83 0.44 -0.94 1.00 0.98 0.93 0.77
Tasa_Intervencion 0.02 0.00 0.02 -0.04 0.92 0.60 -0.86 0.98 1.00 0.85 0.87
Balance_CC 0.00 -0.01 -0.04 -0.12 0.57 0.09 -1.00 0.93 0.85 1.00 0.49
Balance_Fiscal 0.03 0.01 0.08 0.04 1.00 0.91 -0.50 0.77 0.87 0.49 1.00 
corrplot(cor_data, method="circle")

Dado que en este trabajo se pretende realizar la modelación de las variables de costo y gasto de ventas (también considerando la variante en las diferencias de las variables por periodo), la matriz de correlación es un buen indicador para observar relaciones lineales existentes en las variables. En particular, vemos para el caso de la información de gastos, una “fuerte” correlación de las variables de costo, diferencia de costos y diferencia de ventas. Para este caso, es natural encontrar una correlacion positiva con dichas variables, pues es claro que costos y gastos están asociados entre sí. Igualmente, la variable de diferencia de gastos, fue calculada con la variable gastos, por lo tanto tiene sentido encontrar esta correlación. Por otro lado, tanto para la variable de ventas como gastos, no se ven correlaciones fuertes en relación a las variables macro-económicas, lo que nos indica que no existe una relación lineal entre esta variable y las demás. Por lo que una buena alternativa en este trabajo, será considerar transformaciones no lineales de las variables para encontrar una dependencia con la información de los costos de venta. Ahora, mirando las variables de diferencia de gastos y costos, vemos que existen dependencias lineales mayores con las variables macroeconómicas a diferencia de las variables originales, así vemos que la diferencia en costos tiene un grado de asociación con todas las variables macroeconomicas, mientras que los gastos presenta una relación lineal con el balance de cuenta corriente, por lo que se podría pensar que existen relaciones no lineales respecto a los otros indicadores económicos.

Visualizacion de la distribucion de los datos

pairs(data,diag.panel = panel.hist, lower.panel = panel.cor)

Visualizacion de los datos segun curvas de normalidad

### Funcion para encontrar los contornos
c_alpha = function(alpha, sigma, p){
  res = (2*pi)^(-p/2)*(det(sigma))^(-1/2)*exp(-1/2*qchisq(1-alpha, df = p))
  
  return(res)
}

grafica = function(data, name1, name2){

  data_aux = data[c(name1,name2)]
  names(data_aux) = c("y1","y2")
  
  cov_data = cov(data_aux)
  mean_data = mean=colMeans(data_aux)
  
  min_value1 = min(data[name1])
  max_value1 = max(data[name1])
  
  min_value2 = min(data[name2])
  max_value2 = max(data[name2])
  
  n = 100
  
  y1 = seq(min_value1, max_value1, length.out = n)
  y2 = seq(min_value2, max_value2, length.out = n)
  
  grid = expand.grid(y1,y2)
  grid$Z<-apply(grid,1,dmvnorm,mean = mean_data,sigma=cov_data)
  Z<-matrix(grid$Z,nrow=n,ncol=n)
  
  contornos = sapply(c(0.01, 0.05, 0.1), c_alpha, sigma = cov_data, 2)
  
  contour(y1,y2,Z,levels=contornos,labels=c("99%","95%","90%"),
          las=1)
  points(data_aux$y1,data_aux$y2)
  grid()  
  title(main = "Contornos de distribucion normal", xlab = name1, ylab = name2)
}
p1 = grafica(data, "Costo.de.ventas_dif", "PIB")

p2 = grafica(data, "Costo.de.ventas_dif", "TRM")

p3 = grafica(data, "Costo.de.ventas" , "Costo.de.ventas_dif")

p4 = grafica(data, "Gastos.de.ventas" , "Gastos.de.ventas_dif")

p4 = grafica(data, "Costo.de.ventas_dif" , "Gastos.de.ventas_dif")

Creacion de modelos

Para la variable costos de venta, adicional a la transformación que se realizó para obtener la diferencia relativa entre los periodos, en la etapa de modeloamiento, consideraremos para la variable salida la siguiente transformacion:

Ec2 log(1+diff)

Adicionalmente, a los predictores se les realizó la misma transformación.

Agregar un comentario sobre la eleccion de las empresas para train y para test ### Pablo

data_aux = select(all_data,select = -starts_with('NIT'))
data_aux = select(data_aux,select = -starts_with('Year'))

### Transformacion de las variables respuesta

### Recordar para la interpretacion de los modelos:
### regresion log(y) = b*log(x) -> cambio en x, es un aumento en b% en y
### regresion log(y) = b*x -> cambio en x implica un aumento en 100*b puntos en y
### regresion y = b*log(x) -> cambio en x implica un cambio de (b/100)% en y
### regresion y = b*x -> cambio en x implica un aumento de b en y

### Transformacion 1: a nivel logaritmico
data_aux$Costo.de.ventas_dif = log(1 + all_data$Costo.de.ventas_dif)
data_aux$PIB = log(1 + all_data$PIB)
data_aux$TRM = log(1 + all_data$TRM)
data_aux$Desempleo = log(1 + all_data$Desempleo)
data_aux$Balance_CC = log(1 + all_data$Balance_CC)
data_aux$Balance_Fiscal = log(1 + all_data$Balance_Fiscal)
data_aux$Inflacion = log(1 + all_data$Inflacion)
data_aux$Tasa_Intervencion = log(1 + all_data$Tasa_Intervencion)

empresas_train = c(800015615, 800045720,
                   800081030, 800112440,
                   800118660, 800157469,
                   800232356, 800236890,
                   801002644, 805012368,
                   806014553, 830030574,
                   830037495, 860009694,
                   860033653, 860050956,
                   890909034, 900173460,
                   900184722, 900204182,
                   900364670, 900378893)

empresas_test = c(830052054, 860030360,
                   860501682, 890117431,
                   890300012, 890311366,
                   890904459, 890929951,
                   900234565, 900389088,
                   900437650)

train = all_data$NIT%in%empresas_train
test = all_data$NIT%in%empresas_test

NIT_train = NIT[train,]
NIT_test = NIT[test,]

data_train = data_aux[train,]
rownames(data_train) <- NULL
data_train_z = data_train

media_tr = attr(data_train_z, 'scaled:center')
stdev_tr = attr(data_train_z, 'scaled:scale')

data_train_z = as.data.frame(data_train_z)

data_test = data_aux[test,]
rownames(data_test) <- NULL
data_test_z = data_test

data_train_z$NIT = NIT_train
data_test_z$NIT = NIT_test

Definición de las funciones para calcular el R-squared y R-Squared Adjusted

r2_score <- function(x, y) summary(lm(y~x))$r.squared
adj_r2_score <- function(x, y) summary(lm(y~x))$adj.r.squared

Modelo para los Costos

Modelo Lineal General

step.model <- ols_step_both_p(lm('Costo.de.ventas_dif~PIB+TRM+Desempleo+Inflacion+Tasa_Intervencion+Balance_CC+Balance_Fiscal', data = data_train_z), pent = 0.1, prem = 0.1, details = TRUE)
## Stepwise Selection Method   
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. PIB 
## 2. TRM 
## 3. Desempleo 
## 4. Inflacion 
## 5. Tasa_Intervencion 
## 6. Balance_CC 
## 7. Balance_Fiscal 
## 
## We are selecting variables based on p value...
## 
## 
## Stepwise Selection: Step 1 
## 
## - PIB added 
## 
##                           Model Summary                           
## -----------------------------------------------------------------
## R                        0.086       RMSE                  0.457 
## R-Squared                0.007       Coef. Var          -429.959 
## Adj. R-Squared          -0.008       MSE                   0.209 
## Pred R-Squared          -0.045       MAE                   0.240 
## -----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.100         1          0.100    0.478    0.4919 
## Residual       13.367        64          0.209                    
## Total          13.467        65                                   
## ------------------------------------------------------------------
## 
##                                    Parameter Estimates                                    
## -----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig       lower     upper 
## -----------------------------------------------------------------------------------------
## (Intercept)    -0.967         1.246                 -0.776    0.441     -3.456     1.522 
##         PIB    12.709        18.384        0.086     0.691    0.492    -24.017    49.436 
## -----------------------------------------------------------------------------------------
## Warning in if (pvals[minp] <= pent) {: the condition has length > 1 and only the
## first element will be used
## 
## No more variables to be added/removed.
## 
## 
## Final Model Output 
## ------------------
## 
##                           Model Summary                           
## -----------------------------------------------------------------
## R                        0.086       RMSE                  0.457 
## R-Squared                0.007       Coef. Var          -429.959 
## Adj. R-Squared          -0.008       MSE                   0.209 
## Pred R-Squared          -0.045       MAE                   0.240 
## -----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.100         1          0.100    0.478    0.4919 
## Residual       13.367        64          0.209                    
## Total          13.467        65                                   
## ------------------------------------------------------------------
## 
##                                    Parameter Estimates                                    
## -----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig       lower     upper 
## -----------------------------------------------------------------------------------------
## (Intercept)    -0.967         1.246                 -0.776    0.441     -3.456     1.522 
##         PIB    12.709        18.384        0.086     0.691    0.492    -24.017    49.436 
## -----------------------------------------------------------------------------------------
summary(step.model)
##            Length Class      Mode     
## orders      1     -none-     character
## method      1     -none-     character
## steps       1     -none-     numeric  
## predictors  1     -none-     character
## rsquare     1     -none-     numeric  
## aic         1     -none-     numeric  
## sbc         1     -none-     numeric  
## sbic        1     -none-     numeric  
## adjr        1     -none-     numeric  
## rmse        1     -none-     numeric  
## mallows_cp  1     -none-     numeric  
## indvar      7     -none-     character
## betas       2     -none-     numeric  
## lbetas      1     -none-     numeric  
## pvalues     2     -none-     numeric  
## beta_pval   4     data.frame list     
## model      12     lm         list

De los resultados anteriores, vemos que el modelo lineal general con mejor desempeño es aquel que solo considera el PIB además de un intercepto. Pues vemos que en particular para el PIB, su error estandar se reduce considerablemente (un poco más de la mitad), además que vemos mejores valores al R-Squared Adjusted. Por lo tanto, procedemos a mirar el desempeño de un modelo entrenado considerando unicamente el PIB.

mod_lin = lm('Costo.de.ventas_dif~PIB', data = data_train_z)
summary(mod_lin)
## 
## Call:
## lm(formula = "Costo.de.ventas_dif~PIB", data = data_train_z)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.55042 -0.01670  0.08328  0.19822  0.80263 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.9668     1.2460  -0.776    0.441
## PIB          12.7095    18.3843   0.691    0.492
## 
## Residual standard error: 0.457 on 64 degrees of freedom
## Multiple R-squared:  0.007412,   Adjusted R-squared:  -0.008097 
## F-statistic: 0.4779 on 1 and 64 DF,  p-value: 0.4919

De este modelo entrenado, tenemos que un incremento en un 1% del PIB anual del país, corresponde a un aumento del 11.47% anual en la diferencia de costos de venta de las empresas, analogamente al realizar la interpretación sobre la TRM, tenemos que un aumento en un 1% en la TRM anual, implica un aumento del 0.0767% en la diferencia de los costos de venta anuales de las empresas. De acá podemos ver que al utilzar este modelo, el PIB es el que nos presenta un mayor impacto para poder explicar los cambios relativos anuales en los costos de venta de las empresas en el sector de la construcción.

Desempeño del Modelo Entrenado

En el conjunto de entrenamiento:

DEVOVLER LA TRANSFORMACION DEL LOGARITMO

preds = predict(mod_lin)

plot(preds, data_train_z$Costo.de.ventas_dif, xlab = 'Valor predicho', ylab = 'Valor Real', 
     main='Modelo Lineal General')
abline(a=0, b=1, lwd = 2, col = 'red')

En el conjunto de Prueba

DEVOLVER LA TRANSFORMACION DEL LOGARITMO

preds_test = predict(mod_lin, newdata = data_test_z)

r2_model<-r2_score(preds_test, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds_test, data_test_z$Costo.de.ventas_dif)

sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
                sep = " ", collapse = NULL)

plot(preds_test, data_test_z$Costo.de.ventas_dif, xlab = 'Valor predicho', ylab = 'Valor Real',
     main='Modelo Lineal General', sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')

Graficos de Evaluacion de modelos

par(mfrow=c(3,2))
mod_plot = plot(mod_lin, which = c(1:6))

De los gráficos anteriores, en particular analizando el valor de la distancia de Cook para las observaciones, notamos que las observaciones 17, 30, 52 son candidatas a ser registros atípicos en el conjunto de datos, por lo que entrenaremos un nuevo modelo eliminando estos registros. Notar que los outliers corresponden a la empresa con NIT 890909034 para el 2016, y para la empresa 800236890 en los periodos del 2017 y 2018, información que se contrasta con los outliers obtenidos en el boxplot, lo cual están incluidos tanto mirando la variable de diferencia costos de venta.En vista de su comportamiento raro en magnitud según el boxplot, no es de extrañarse que estos registros hayan sido aquellos con la distancia de cook mas grande. Procederemos a entrenar un modelo removiendo estos registros

### Vector para eliminar registros atipicos según la información obtenida por la distancia de cook
outliers = c(17, 30, 52)
data_train_z_noOut = data_train_z[-outliers,]
h_ii<-hatvalues(mod_lin)
plot(hatvalues(mod_lin),las=1,xlab="i",ylab="hii",main="Influencia (h_ii)",type="h")

### Inflacion de la varianza

#vif(mod_lin)

Entrenamiento del modelo sin Outliers

mod_lin02 = lm('Costo.de.ventas_dif~PIB', data = data_train_z_noOut)
summary(mod_lin02)
## 
## Call:
## lm(formula = "Costo.de.ventas_dif~PIB", data = data_train_z_noOut)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83980 -0.08663  0.00809  0.11699  0.69910 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -1.3702     0.6017  -2.277   0.0263 *
## PIB          19.8623     8.8777   2.237   0.0289 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2156 on 61 degrees of freedom
## Multiple R-squared:  0.07584,    Adjusted R-squared:  0.06069 
## F-statistic: 5.006 on 1 and 61 DF,  p-value: 0.02893

Después de remover los outliers, tenemos que un incremento en un 1% del PIB anual del país, corresponde a un aumento del 32.81% anual en la diferencia de costos de venta de las empresas (note que ahora el PIB tiene un impacto mayor que en el modelo anterior, pues se nota la diferencia en su coeficiente, un aumento casi de 3 veces el valor del modelo anterior), analogamente al realizar la interpretación sobre la TRM, tenemos que un aumento en un 1% en la TRM anual, implica un aumento del 0.07627% en la diferencia de los costos de venta anuales de las empresas (lo cual mantiene un comportamiento similar al indicado por el modelo anterior). De acá podemos ver igualmente que al utilzar este modelo, el PIB es el que nos presenta un mayor impacto para poder explicar los cambios relativos anuales en los costos de venta de las empresas en el sector de la construcción.

Desempeno del Modelo Entrenado sin Outliers

En el conjunto de entrenamiento:

preds = predict(mod_lin02)

plot(preds, data_train_z_noOut$Costo.de.ventas_dif, xlab = 'Valor predicho', ylab = 'Valor Real')
abline(a=0, b=1, lwd = 2, col = 'red')

En el conjunto de Prueba

preds_test = predict(mod_lin02, newdata = data_test_z)

r2_model<-r2_score(preds_test, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds_test, data_test_z$Costo.de.ventas_dif)

sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
                sep = " ", collapse = NULL)

plot(preds_test, data_test_z$Costo.de.ventas_dif, xlab = 'Valor predicho', ylab = 'Valor Real',
     main='Modelo Lineal General', sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')

Graficos de Evaluacion de modelos entrenando sin Outliers

par(mfrow=c(3,2))
mod_plot = plot(mod_lin02, which = c(1:6))

Al observar el desempeño del modelo sin considerar dichos registros raros, no se tiene una mejora significativa en el modelo.

Modelos Lineales Generalizados (GLM)

Dentro de los modelos que fueron expuestos a lo largo del curso, se encuentra la familia de los modelos lineales generalizados, los cual nos permiten la creación de modelos teniendo en cuenta diferentes supuestos que se hacen sobre la variable respuesta que estamos considerando. En la elaboración de este trabajo, se está considerando realiza un modelo de predicción sobre transformaciones en la diferencia relativa de los costos de ventas que tienen las empresas, por lo que el dominio de la variable son los números reales (ya que puede tomar valores continuos tanto positivos como negativos, los cuales son considerados en las diferentes transformaciones que se plantean). Justo por la caracteristica de la variable respuesta que tenemos, no se considera pertinente realizar evaluacion de los modelos, esta conclusión está apoyada del análisis que se realiza de los diferentes modelos generalizados que podriamos considerar:

Por el análisis realizado, no se considera pertinente utilizar alguno de los otros modelos lineales generalizados para el modela miento de la variable transformada de diferencia relativa de costos de venta

Modelo de Efectos Mixtos

mod_me = lmer('Costo.de.ventas_dif~Desempleo+PIB+(0+Desempleo|NIT)+(PIB|NIT)', data = data_train_z)
summary(mod_me)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Costo.de.ventas_dif ~ Desempleo + PIB + (0 + Desempleo | NIT) +  
##     (PIB | NIT)
##    Data: data_train_z
## 
## REML criterion at convergence: 68.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2577 -0.0832  0.1859  0.3549  2.5374 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  NIT      Desempleo   3.576169 1.89108       
##  NIT.1    (Intercept) 0.004193 0.06475       
##           PIB         0.566302 0.75253  -1.00
##  Residual             0.182955 0.42773       
## Number of obs: 66, groups:  NIT, 22
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -0.8466     3.0990  -0.273
## Desempleo    -1.2765    30.4740  -0.042
## PIB          12.6315    17.3072   0.730
## 
## Correlation of Fixed Effects:
##           (Intr) Desmpl
## Desempleo -0.926       
## PIB       -0.473  0.107

Ahora considerando el modelo de efectos mixtos, observando primero la componente de los efectos fijos, vemos que un incremento en un 1% del PIB anual del pais, refleja un aumento en un 12.63% en la diferencia de los costos de venta anuales de las empresas, mientras que al ver un cambio de un 1% en el desempleo, este disminuye la diferencia de los costos en un 1.2765%, ahora bien, considerando los efectos aleatorios, los cuales tendrán un efecto sobre la predicción en aquellas empresas que fueron consideradas en el entrenamiento (mirando al NIT), vemos que adicional a la componente fija, el cambio de 1% del PIB anual, que antes tenia un efecto en el aumento del 12.63% sobre la diferencia de los costos de venta, ahora se aumenta en un 0.56% adicional, es decir, para aquellas empresas que haya visto el modelo previamente, el aumento anual en los costos de venta corresponderá a un 13.198%. Ahora realizando el análisis de los efectos aleatorios para el desempleo, vemos que el desempleo tiene un incremento del 3.576% en la diferencia de los costos de venta anual por un aumento del 1% en el desempleo, por lo que para aquellas empresas consideradas al momento de entrenar el modelo, presentarán un aumento del 2.299669% en la diferencia de costos anuales por un cambio de un 1% en el desempleo, por lo que estas empresas pueden tener caracteristicas que evitan que el desempleo no tenga un imacto muy grande en sus costos.

par(mfrow=c(3,2))
plot(mod_me, which = c(1:6))

DEVOLVER TRANSFORMACIONES DEL LOGARITMO

preds = predict(mod_me)
plot(preds,data_train_z$Costo.de.ventas_dif, main='Modelo de Efectos Mixtos')
abline(a=0, b=1, lwd = 2, col = 'red')

DEVOLVER TRANSFORMACIONES DEL LOGARITMO

preds = predict(mod_me, newdata = data_test_z,allow.new.levels=TRUE)

r2_model<-r2_score(preds, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds, data_test_z$Costo.de.ventas_dif)

sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
                sep = " ", collapse = NULL)

plot(preds,data_test_z$Costo.de.ventas_dif, main='Modelo de Efectos Mixtos', sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')

Desafortunadamente en el caso de prueba no se logra distinguir la tendencia natural en el comportamiento de la diferencia de los costos de venta y la escala se encuentra muy distante entre los valores predichos y los reales, como se espera desde el gráfico anterior. En esta evaluación los valores predichos se ponderan para generar una tendencia casi constante y la razonabilidad una vez más se pone en duda, aunque no por completo, dado que éstos valores efectivamente pueden ser alcanzados con condiciones ligeramente diferentes.

Superficie poblacional del modelo

Desempleo = seq(min(data_train_z$Desempleo),max(data_train_z$Desempleo), length.out = 50)
PIB = seq(min(data_train_z$PIB),max(data_train_z$PIB), length.out = 50)

data_surface = expand.grid(Desempleo,PIB)
names(data_surface) = c('Desempleo','PIB')
data_surface$NIT = 1000

z = predict(mod_me, newdata = data_surface,allow.new.levels=TRUE)
z = matrix(z, nrow = 50, ncol = 50)

### Toca rotar la figura para poder ver el plano
fig <- plot_ly(x = PIB, y = Desempleo, z = z) %>% add_surface()
fig

Superficie del modelo para un individuo conocido

Desempleo = seq(min(data_train_z$Desempleo),max(data_train_z$Desempleo), length.out = 50)
PIB = seq(min(data_train_z$PIB),max(data_train_z$PIB), length.out = 50)

data_surface = expand.grid(Desempleo,PIB)
names(data_surface) = c('Desempleo','PIB')
data_surface$NIT = 800015615

z2 = predict(mod_me, newdata = data_surface)
z2 = matrix(z2, nrow = 50, ncol = 50)

### Toca rotar la figura para poder ver el plano
fig %>% add_surface(z = ~z2, opacity = 0.98,colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(120,0,64)"))) %>% layout(title="Poblacional (azul y verde) vs Un Individuo (Rosa)")
#fig <- plot_ly(x = PIB, y = Desempleo, z = z) %>% add_surface()
#fig

AGREGAR SIN OUTLIERS EN MODELO MIXTO